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MS.  ABSTRACT 


— Two-dimensional  numerical  techniques  have  been  applied  to  obtain  solutions  of 
problems  of  stress  wave  interactions  with  cracks  in  rock  media.  The  SHEP  code,  a 
finite-difference  Iagrangian  program  that  incorporates  a  comprehensive  elastic-plastic- 
hydrodynamic  behavioral  model,  was  utilized  as  the  basic  numerical  program  in  this 
study.  Significant  new  capabilities  of  the  code  developed  during  this  program  were  (a) 
the  treatment  of  antiplane  shear,  or  out-of-plane  displacement  (where  the  motion  is  de¬ 
pendent  only  on  the  in-plane  coordinates)  and  (b)  the  incorporation  of  a  dynamic  brittle 
fracture  criterion, (known  as  the  Griffith  criterion). 

CyO 

In-plane  problems  considered  were  the  propagation  of  a  stress  wave  through  a 
SG!w"\n^n^e  sPa<^e  containing  (a)  a  semi-infinite  crack,  (b)  a  finite  stationary  crack,  1 
and  (c)  a  finite  running  crack.  The  crack  was  oriented  at  a  30°  angle  with  the  wave  from 
in  these  problems.  The  computations  demonstrated  that  it  is  possible  to  observe  the 
diffraction  of  the  shear  wave  fronts  and  to  follow  the  flow  of  energy  through  the  crack. 

An  ideal  but  relatively  simple  problem  was  chosen  as  a  test  vehicle  for  the  ouc-of- 
plane  program,  namely  the  antiplane  shear  loading  of  a  crack.  In  the  case  of  a  station¬ 
ary  crack,  the  code  solution  was  shown  to  give  excellent  agreement  in  all  respects  with 
the  analytical  solution.  Using  the  Griffith  criterion,  moderately  good  agreement  was 
obtained  for  the  case  of  the  running  crack. 
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ABSTRACT 


Two-dimensional  numerical  techniques  have  been  applied  to  obtain 
solutions  of  problems  of  stress  wave  interactions  with  cracks  in  rock  media. 
The  SHEP  code,  a  finite-difference  Lagrangian  program  that  incorporates  a 
comprehensive  elastic -plastic-hydrodynamic  behavioral  model,  was  utilized 
as  the  basic  numerical  program  in  this  study.  Significant  new  capabilities 
of  the  code  developed  during  this  program  were  (a)  the  treatment  of  antiplane 
shear,  or  out-of-plane  displacement  (where  the  motion  is  dependent  only  on 
the  in-plane  coordinates)  and  (b)  the  incorporation  of  a  dynamic  brittle  frac¬ 
ture  criterion  (known  as  the  Griffith  criterion). 

In-plane  problems  considered  were  the  propagation  of  a  stress  wave 
through  a  semi-infinite  space  containing  (a)  a  semi-infinite  crack,  (b)  a 
finite  stationary  crack,  and  (c)  a  finite  running  crack.  The  crack  was 
oriented  at  a  30°  angle  with  the  wave  front  in  these  problems.  The  compu¬ 
tations  demonstrated  that  it  is  possible  to  observe  the  diffraction  of  the 
shear  wave  fronts  and  to  follow  the  flow  of  energy  through  the  crack. 

An  ideal  but  relatively  simple  problem  was  chosen  as  a  test  vehicle 
for  the  out-of-plane  program,  namely  the  antiplane  shear  loading  of  a  crack. 
In  the  case  of  a  stationary  crack,  the  code  solution  was  shown  to  give 
excellent  agreement  in  all  respects  with  the  analytical  solution.  Using  the 
Griffith  criterion,  moderately  good  agreement  was  obtained  for  the  case  of 
the  running  crack. 

In  running  the  SHEP  program,  it  is  necessary  to  dampen  oscillations 
at  the  wave  front  by  the  incorporation  of  a  fictive  viscosity.  In  the  case  of 
pressure,  a  quadratic  viscosity  component  is  added.  In  the  case  of  shear, 
a  linear  term  is  normally  added.  A  linear  viscosity  tends  to  increasingly 
spread  the  wave  front,  however.  A  study  of  the  origins  of  this  viscosity 
term  has  led  to  the  development  of  an  alternate  formulation  for  the  viscosity 
term  which  is  cubic  and  applies  both  to  pressure  and  shear.  This  alternate 
formulation  will  eventually  replace  the  so-called  anisotropic  or  Navier 
Stokes  linear  viscosity. 
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INTRODUCTION 


The  objective  of  this  program  was  to  apply  two-dimensional  numer¬ 
ical  techniques  to  the  solution  of  problems  of  stress  wave  interactions  with 
individual,  plane  cracks  in  rock  media.  Computer  codes  suitable  for  solving 
a  wide  range  of  fluid  and  solid  mechanics  problems  have  been  available  cr 
under  development  for  several  years  and  efforts  to  extend  these  techniques 
to  the  quantitative  analysis  of  wave  interactions  with  cracks  and  of  crack 
motion  are  now  feasible.  Shock  Hydrodynamics  two-dimensional  SHEP  (Shock 
Hydrodynamic  Elastic  Plastic)  code  was  utilized  for  this  purpose  in  this 
program.  Analytical  studies  were  also  conducted  in  conjunction  with  the 
numerical  work. 

2.  '  SUMMARY 

2.1  PROBLEM  AREA 

Excavation  processes  in  rock  media  typically  involve  the  action 
of  strong  dynamic  stresses  introduced  either  from  explosive,  mechanical,  or 
other  impulsive  loading  sources.  The  propagation  of  stress  waves  in  homo¬ 
geneous,  isotropic  media  is  reasonably  well  understood.  In-situ  rock  media, 
however,  typically  contain  large  scale  discontinuities  in  the  form  of  cracks, 
joints,  and  faults.  Interactions  of  stress  waves  with  these  discontinuities 
can  cause  slippage  along  the  cracks,  or  extension  (propagation)  of  the 
cracks,  or  separation.  The  interactions  can  also  alter  the  characteristics 
of  the  stress  wave  transmitted  across  the  discontinuity. 

This  study  is  concerned  with  an  analysis  of  the  detailed  mechan¬ 
isms  involved  when  strong  stress  waves  interact  with  the  crack  surfaces  in 
jointed  rock  media.  Of  particular  interest  are  cracks  which  are  obliquely 
oriented  relative  to  the  wave  front.  The  understanding  thus  obtained  can 
contribute  to  the  advancement  of  knowledge  of  excavation  processes  in 
various  media,  and  how  to  control  and/or  improve  such  processes. 
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The  technical  approach  used  to  study  the  details  of  stress  wave-  i 
crack  Interactions  was  based  on  two-dimensional  numerical  analyses  of  the 
dynamic  phenomena  occurring  under  various  conditions  of  stress  wave  loading 
and  crack  orientation  relative  to  the  wavd  front.  The  computer  program  used 
to  obtain  the  numerical  solutions  was  the  SHEP  code.  SHEP  is  a  finite- 
difference  Lagrangian  program  employing  a  comprehensive  hydrodynamic- 
e la stlc -plastic  behavioral  model.  SHEP  has  been  under  intensive  use  and 
development  for  the  past  six  years  and  has  been  applied  to  a  broad  spectrum 
of  wave  propagation  problems.  , 

A  major  difficulty  in  examining  wave  interactions  with  discontinui¬ 
ties  such  as  cracks  or  fracture  surfaces. arises  due  to  the  constraint  of  the 
continuum  model  which  is  normally  assumed  in  numerical  analyses  of  wave 
propagation.  Special  routines  in  the  SHEP  code  alleviate  this  difficulty  by 
permitting  crack  surfaces  to  be  explicitly  defined  in  the  computational  grid,. 
Thus  the  grid  is  not  coupled  across  the  crack,  and  slippage  and/or  separa¬ 
tion  can  occur.  1 

2.2  PLAN  OF  RESEARCH  •  , 

The  contract  work  was  divided  into  the  following  two  ^sks, 

corresponding  to  analyses  of  stress  wave  interactions  with  1  , 

\  ,  ; 

a)  single,  semi-infinite  cracks  and 

b)  single,  finite -length  cracks.  , 

,  } 

2.2.1  Task  1  -  Interaction  of1  Stress  Waves  with  Single. 

Semi -Infinite  Cracks 

This  task  was  concerned  with  the  analysis  pertaining  to  cracks 
which  are  semi -infinite  in  extent,  or  which  intersect  with  the  ground  sur- 

I 

face.  The  work  initially  consisted  of  the  selection  and  determination  of 
the  problem  specifications  (such  as  media  properties,  loading  wave  char¬ 
acteristics,  and  crack  orientation  and  condition).  SHEP  code  solutions  of 
the  problems  defined  were  then  set  up  and  run.  These  solutions  provide 
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complete  quantitative  data  for  all  the  state  and  motion  variables  of  interest 
throughout  the  computing  field  at  regular  intervals  of  time.  In  addition, 
spatial  plots  of  the  principai stress,  particle  velocity,  and  displacement 
fields  were  obtained  a^  selected  times  during  the  event.  Following  comple¬ 
tion  of  the  solutions,  analysis  and  interpretation  of  the  results  were  per¬ 
formed  to  characterize  the  transmitted  stress  waves. 

r  > 

i  i  1 

| 

i  2.2.2  Task  2  -  Interaction  of  Stress  Waves  with  Single, 
i  Finite  Cracks 

rr> . - 

This  task  was  concerned  with  cracks  of  finite  length,  so  that  the 
interaction  of  a  stress  wave  with  a  crack  tip  could  be  examined.  The  prob¬ 
lems  wfere  selected  so  that  comparisons  between  the  cases  of  semi -infinite 
(Task  1)  and  finite  (Task  2)  cracks  could  be  made. 

1  Also  included  in  this  task  was  the  development  of  a  dynamic 
Griffith  fracture  criterion,  to  enable  the  study  of  crack  propagation.  In 
conjunction  with' this,  the  SHEP  code  was  generalized  to  accomodate  anti¬ 
plane  shear,  or  out-of-plane  displacement. 

T 6  verify  the  suitability  and  accuracy  of  the  code,  and  modifica¬ 
tions  made  thereto  during  the  course  of  the  program,  for  analyses  of  wave/ 
crack  interactions  ,  analytic  solutions  of  selected  test  problems  were  obtained 
and  used  to  check  out  corresponding  numerical  solutions  obtained  with  the 
i  code.  1  , 

j 

2.3  major  Accomplishments 

I  ! 

Pribary  program  accomplishments  included 

a)  The  completion  of  SHEP  code  solutions  of  three 
1  in-plane  stress  wave/crack  interaction  problems. 

j  b)  The  codification  of  the  equations  of  anti-plane 
1  shear  (z -independent)  to  provide  a  new  capability 
for  SHEP  code  analyses. 

■  i  i  1 

,  c)  The  formulation  of  a  dynapic  Griffith  criterion 

1  for  analyses  of  crack  propagation  with  the  SHEP  code. 
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These  items  are  summarized  in  the  following  sub-sections. 

2.3.1  In-Plane  Numerical  Solutions 

Three  problems  of  wave  interactions  with  cracks  were  selected  for 
analysis  by  means  of  code  solution,  as  depicted  in  Figure  1.  These  problems 
were  chosen  to  demonstrate  the  utility  of  numerical  techniques  for  obtaining 
detailed  information  on  the  response  of  cracked  media  subjected  to  Impulsive 
loads,  in  general,  and,  in  particular,  for  assessments  of  the  wave  inter¬ 
actions  in  the  vicinity  of  a  crack.  It  is  noted  that  the  explicit  definition  of 
a  crack  surface  such  as  specified  in  these  problems  is  not  generally  amen¬ 
able  to  treatment  through  conventional  code  techniques,  which  normally 
assume  a  continuous  material  model. 

The  rock  medium  selected  for  these  problems  was  granite.  The 
material  properties  assumed  for  the  granite  are  described  in  Section  3.2  of 
this  report.  The  problems  were  run  in  plane  geometry,  assuming  plane  strain; 
the  variables  are  thus  independent  of  the  z-coordinate  (perpendicular  to  the 
cross-section  shown). 

The  first  problem  selected  for  analysis  (Case  1) -consisted  of  inter¬ 
actions  of  a  stress  wave  with  a  crack  oriented  at  a  30°  angle  with  the  wave 
front.  The  stress  wave  was  generated  by  uniformly  loading  the  left  face  of 
the  granite  block  with  a  pressure  pulse.  A  triangular  pressure  pulse  of  5 
kllobars  peak  magnitude  and  0.2  millisecond  duration  was  used  for  this 
problem,  as  sketched  below: 


P  (kb) 
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CASE  1.  INTERACTION  OF  STRESS  WAVE  WITH  SINGLE,  SEMI -INFINITE  CRACK 


CASE  2.  INTERACTION  OF  STRESS  WAVE  WITH  SINGLE,  FINITE-LENGTH  CRACK 


CASE  3.  INTERACTION  OF  STRESS  WAVE  WITH  SINGLE,  FINITE-LENGTH 
CRACK,  WITH  CRACK  GROWTH 


Figure  1.  Specifications  of  Problems  for  SHEP  Code  Solutions 


The  crack  in  this  problem  was  characterized  by  a  free  slip  condi¬ 
tion  and  zero  width.  Opening  of  the  crack  was  allowed  to  occur  if  stress 
components  normal  to  the  crack  went  into  tension. 

The  second  problem  considered  (Case  2)  involved  the  interaction 
of  a  stress  wave  with  a  finite-length  crack.  The  angle  of  orientation  of  the 
crack  and  the  applied  pressure  loading  were  the  same  as  in  the  first  problem. 
The  crack  extended  from  the  loading  surface  (lower  left)  to  the  crack  tip, 
situated  on  the  horizontal  mid-plane  of  the  block. 

Case  3  was  the  same  as  the  second  problem,  except  that  crack 
growth  was  permitted.  This  case  demonstrates  the  provisions  in  the  code 
which  can  be  used  to  model  crack  propagation.  In  this  case,  dynamic  de¬ 
coupling  of  lattice  points  in  the  computing  mesh  occurs  when  a  specified 
criterion  is  satisfied. 

SHEP  code  solutions  of  these  three  problems  were  successfully 
completed.  Plots  depicting  the  particle  velocity  field  and  the  principal 
stress  field  occurring  in  the  test  block  were  obtained  for  several  times  during 
the  interactions.  In  addition,  time  histories  of  pertinent  parameters  at 
several  stations  in  the  field  were  recorded.  These  results  are  discussed  in 
detail  in  Section  3  of  the  report.  Some  representative  results  of  these  code 
solutions  are  shown  here,  in  Figures  2  to  7. 

For  Case  1,  the  principal  stress  field,  for  a  time  of  .3  msec,  and 
the  particle  velocity  field,  for  a  time  of  .5  msec,  are  shown  in  Figures  2 
and  3.  As  the  wave  encounters  the  crack,  the  principal  stress  vectors  may 
be  seen  (Figure  2)  to  rotate  into  a  direction  transverse  to  the  crack  surface, 
reflecting  the  fact  that  the  crack  surface  can  not  bear  shear  stress.  This 
interaction  produces  a  dilatational  wave  and  trailing  shear  wave  which 
propagate  across  the  block,  as  indicated  in  the  velocity  field  plot  (Figure  3). 
The  peak  stress  in  the  transmitted  wave  was  reduced  by  about  25%  from  that 
in  the  incident  wave. 


6 


8H0CN  HVOMOVKftMCS.  INC.  SHCf  coor 
RUN  NO.  7119-1.  MflVC  INTCIRCT1IN  H17N  9O-0CO  CMC* 


CYCUC  51 


T  =  502.5B61  USfC 


Figure  3.  Particle  Velocity  Field,  Case  1 
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Figure  5.  Particle  Velocity  Field,  Case  3 


Figure  6.  Particle  Velocity  Field,  Case  3. 


An  example  of  the  results  of  the  Case  2  solution  is  shown  in 
Figure  4,  which  depicts  the  particle  velocity  field  occurring  at  a  time  of 
.6  msec.  In  this  interaction,  the  wave  system  is  divided  approximately  in 
half,  the  part  above  the  crack  tip  appearing  as  a  simple  plane  wave,  and 
that  below  as  a  dilate tional  wave  and  trailing  shear  wave,  as  in  the  pre¬ 
vious  solution.  Starting  from  the  crack  tip,  a  disturbance  propagates  into 
the  plane  wave  region  above  and  the  "cracked"  region  below,  altering  both 
flow  fields  and  creating  an  expanding  region  of  transition  between  them. 

In  Case  3,  where  crack  growth  was  allowed,  the  criterion  used 
for  decoupling  of  points  beyond  the  crack  tip  was  that  each  of  the  cells 
surrounding  a  lattice  point  must  have  failed,  i.e. ,  at  some  time  reached  a 
state  on  the  granite  failure  surface.  This  was  a  conservative  criterion,  since 
the  failure  surface  generally  represents  states  where  virtually  complete  frac¬ 
ture  occurs.  Use  of  more  sensitive  criteria  can  be  employed,  and  develop¬ 
ment  of  a  dynamic  Griffith  criterion  was  subsequently  undertaken,  as 
described  below.  The  crack  growth  occurring  in  this  solution  is  indicated 
in  Figures  5,  6,  and  7,  which  are  plots  of  the  particle  velocity  field  for 
times  of  .5,  .6,  and  .7  msec.  The  extent  of  the  crack  in  each  plot  is  indi¬ 
cated  by  the  circled  lattice  points. 

2.3.2  Dynamic  Griffith  Criterion 

In  connection  with  efforts  made  under  the  program  to  enable  the 
study  of  propagation  of  cracks  under  stress  wave  loading,  the  incorporation 
of  equations  into  the  SHE P  code  which  govern  the  rate  of  propagation  of  a 
brittle  crack  surface  in  an  elastic  material  was  undertaken.  These  equations 
are  known  as  the  Griffith  criterion  and  they  provide  a  relation  between  power 
input  to  the  body  and  the  rate  of  uptake  of  this  power  by  strain  energy, 
kinetic  energy,  and  new  surface  energy.  The  concept  of  surface  energy  is 
the  feature  that  was  introduced  by  Griffith  in  the  early  1900's,  and  it 
requires  the  determination  of  an  additional  material  parameter,  namely  the 
surface  energy  per  unit  area.  An  algorithm  appropriate  for  the  SHEP  code 
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was  programmed  and  tested  on  a  model  problem.  The  results  indicated  that 
the  effective  crack  propagation  was  slower  than  predicted  theoretically.  It 
is  expected  that  Improved  results  would  be  obtained  with  an  alternate  form 
of  the  flctive  viscosity. 

2.3.3  Analytical  Comparison  Problems 

To  verify  the  SHIP  code  solutions  and  the  formulation  changes 
made,  comparisons  of  numerical  results  with  analytical  solutions  of  model 
problems  were  made. 

As  part  of  this  effort,  the  capability  was  added  to  the  code  for 
the  treatment  of  anti-plane  shear,  or  out-of-plane  displacements,  with  the 
restriction  that  the  motions  are  Independent  of  the  z-coordlnate,  so  as  to 
retain  the  two-dimensional  character  of  the  code.  This  was  done  primarily 
since  the  only  elasto-dynamlc  solutions  currently  available  for  an  accelera¬ 
ting  crack  are  those  for  the  case  of  anti-plane  shear,  although  it  also  repre¬ 
sents  a  useful  tool  in  numerical  analysis  which  has  heretofore  been 
unavailable. 

A  model  problem  of  simple,  shear  motion  of  a  slab  was  first 
solved  with  the  modified  code.  The  results  of  the  code  solution  showed 
excellent  agreement  with  the  analytical  solution  for  this  case.  A  full,  two- 
dimensional  problem  involving  the  interaction  of  an  anti -plane  shear  wave 
with  a  stationary  crack  was  then  set  up  and  run.  Excellent  agreement  with 
the  analytical  solution  was  also  achieved  for  this  case. 

3.  NUMERICAL  SOLUTIONS  OF  IN-PLANE  PROBLEMS 

As  noted  above,  numerical  solutions  of  three  problems  Involving 
the  interaction  of  stress  waves  with  cracks  were  performed.  The  specifica¬ 
tions  cf  these  problems  were  given  in  Section  2.3.1. 
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3.1 


COMPUTATIONAL  METHOD 


3.1.1  Physical  Model 

The  computer  program  used  in  this  study  was  the  two-dimensional 
SHEP  code,  which  solves  the  equations  of  motions  for  elastic-plastic  bodies 
by  means  of  a  finite-difference  Lagrangian-cell  technique.  SHEP  has  been 
under  intensive  use  and  development  for  the  past  six  years  and  has  been 
previously  documented1  and  distributed  to  interested  parties.  The  mathema¬ 
tical  formulation  is  basically  the  same  as  that  described  by  Wilkins9.  To 
delineate  the  boundary  between  elastic  and  plastic  deformations,  various 
yield  criteria  may  be  used,  such  as  von  Mlses.  Mohr-Coulomb,  or  arbitrary 
functions.  Within  the  chosen  yield  surface ,  the  deformations  are  con¬ 
sidered  to  be  elastic,  i.e. ,  when 

J*T2  <  Y  (1) 

where  is  the  second  invariant  of  the  deviatoric  stress  tensor  and  Y  is  the 
yield  strength.  Excursions  on  the  yield  surface  can  be  made  in  accordance 
with  either  the  Prandtl-Reuss  or  plastic  potential  flow  rules. 

To  model  a  crack  surface,  SHEP  contains  provisions  for  inserting 
surfaces  of  discontinuity,  which  consist  of  grid  lines  having  a  dual  set  of 
lattice  points.  These  surfaces  are  discussed  in  the  following  section. 

3.1.2  Surfaces  of  Discontinuity 

In  a  normal  Lagrangian  computational  grid,  material  elements  on 
either  side  of  an  inter  face  at  any  point  are  coupled  to  each  other  for  the 
entire  problem;  they  are,  i.e. ,  locked  or  welded  together  along  the  line 
segment  connecting  any  two  lattice  points  along  the  interface.  At  any  inter¬ 
face,  which  may  represent  a  crack  within  a  material  or  the  boundary  between 
two  different  materials,  there  are,  however,  in  general,  special  boundary 
conditions  which  apply,  and,  in  addition,  there  is  the  possibility  of  forces 
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which  may  be  set  up  that  tend  to  cause  the  materials  to  slip  past  each 
other  or  to  separate.  A  gas  flowing  past  a  metal  surface  is  an  ekample  of 
such  a  case.  The  onset  of  material  fracture  during  a  problem  also  gives  rise 
to  the  requirement  for  treating  the  decoupling  or  uncoupling  of  elements 
which  are,  in  this  case,  within  an  originally  competent  material.  For  appli¬ 
cation  to  problems  in  fracture  mechanics,'  such  as  !ln  this  program,  ttye  latter 
requirement  is  particularly  Important.  ‘ 

t 

1  ,  I 

A  formulation  of  sliding  interfaces  for  Lagranglan  codes,  as  re¬ 
ported  by  Wilkins3,  provided  a  capability  for  the  numerical  treatment  of 

problems  involving  sliding  of  two  materials  along  an  Interface,  ^hls  formu- 

1  , 

latlon  served  as  the  basis  for  develop.,  ent  of  the  surface  of  discontinuity 

*  i 

capability  currently  available  in  the  SHEf  code.  i 

The  basic  features  of  the  surface  of  discontinuity  formulation  are 
illustrated  in  Figure  8.  The  grid  line  corresponding  to  the  surface  of  dis-, 
continuity  is  known  in  common  parlance  as  a  slide  line.  At  the  start  bf  a  , 
problem,  the  lattice  points  along  the  slide  line  mhy  be  individually  desig- 

i  I 

nated  as  decoupled  points,  corresponding  to  their  lying  on  an  interface,  or 

as  coupled  points,  in  which  case  their  behavior  is  the  same  as  in  an  ordinary 

mesh.  For  decoupled  points,  special  sets  of  governing  equations  are  used 

to  individually  determine  the  motion  of  the  point  pairs,  to  reflect  the  fact 

that  there  is  an  interface,  such  that,  e.g. ,  shear  stress  cannot  be  suppprted. 

If  forces  are  present  which  tend  to  cause  slippage,  the  decoupled  points  will 

»  |  • 
thus  disengage  and  move  separately  along  the  slide  line.  , 

Additionally,  the  development  of  tensile  stresses  normal  to  an 
Interface  will  tend  to  cause  material  separation  and  formation  of  voids.  Pro¬ 
visions  have  been  made  in  the  code  to  treat  this  phenomenon,  also. 

I 

The  void  opening  test  is  made  by  computing  the  stress  pormal  to 
the  interface  at  a  decoupled  slide  point  and  comparing  this  value  with  a 
selected  critical  value  of  stress  required  for  uncoupling.  If  the  computed 
stress  is  greater  than  the  critical  value  (in  tension),  then  that  point  is 
designated  as  a  free  surface  point.  The  newly  formed  free  point  is  then  . 
moved  in  accordance  with  the  regular  equations  of  motion  for  a  point  on  a 
free  surface.  1  i  1 
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a.  Normal  Coupled . 
Lagrangian  jGrid 


b .  Decoupled  Grid  With  Relative 
Motion  (Slippage)  Along  Material 
Interface,  Fracture  Surface,  or 
Plane  of  Weakness 


.  "Crack  Propagation"  With  d.  "Crack  Propagation"  With 

Subsequent  Slippage  Along  Subsequent  Separation  and 

Fracture  Surface  Rejoining  Along  Fracture  Surface 

,  1  ' 


Figure  8.  Schematic  Illustrations  of  the  Treatment  of  Slippage, 
Fracture Void  Opening,  and  Void  Closing  Along 
Surfaces  of  Discontinuity  jn  the  SHEP  Code. 


i 


i 


! 


17 


For  the  problems  performed  In  this  study,  the  critical  value  for 
uncoupling  was  set  to  zero,  such  that  any  tensile  stress  would  tend  to  cause 
separation.  In  other  applications,  e.g.,  for  the  Interfaces  In  laminated 
materials,  this  value  could  be  set  equal  to  the  bond  strength. 

Void  closure  may  also  occur  and  is  treated  in  the  code  by  appro¬ 
priate  tests  to  determine  if  the  materials  have  come  in  contact.  If  so,  the 
equations  of  the  surface  of  discontinuity  are  restored,  and  the  materials  may 
subsequently  slip  or  re-open,  as  before. 

Lattice  points  along  the  slide  line  which  are  initially  designated 
as  coupled  points  may  dynamically  decouple,  individually,  during  the  course 
of  a  problem,  if  a  selected  criterion  is  met.  Various  decoupling,  or  fracture, 
criteria  may  be  used.  Once  the  lattice  points  are  decoupled,  the  equations 
of  the  surface  of  discontinuity  are  invoked.  A  mechanism  is  thus  provided 
which  can  be  used  to  model  crack  propagation  within  a  material. 

Additional  information  on  the  code  operations  and  the  mathema¬ 
tical  basis  of  sliding  interfaces  has  been  previously  reported.3”4 

3.2  MATERIAL  PROPERTIES 

The  rock  medium  in  these  problems  was  granite.  The  properties 
selected  for  the  granite  were: 

3 

Density:  pQ  *  2.69  gm/cm 

Dilatational  Velocity:  v,  =  .579  cm/(isec 

ao 

Shear  Velocity:  v_  =  .  330  cm/|isec 

so 

These  values  imply  the  following  other  properties: 

Bulk  Modulus:  KQ  =  .512  Mb 

Shear  Modulus:  G0  =  .293  Mb 

Poisson's  Ratio:  vQ  =  .26 
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The  subscript  o  in  the  above  indicates  that  these  are  normal,  pre-shocked 
values.  These  values  were  selected  from  previous  studies  involving  granite 
media.6*6 . 

The  equation  cf  state  of  granite,  suitable  for  the  low-pressure 
regime  applicable  in  these  problems,  was  formulated  as  follows. 

P  =  A|a  +  Bn2  6  +  Gpe  P<  .04  Mb  (2) 

where 

6  =  1  for  n  >  0 
6  =  0  for  n  £  0 

The  symbols  are  defined  as 

e  =  specific  internal  energy 
P  =  pressure 

r\  =  p/PQ  =  relative  density 
|i  =  r\  -  1  =  compression 
p  =  density 

The  values  of  the  coefficients  are: 

A  =  .512  Mb 

B  =  1.49  Mb 
G  =  2.1 

No  hydrostatic  tension  was  permitted,  i.e.,  P|njn  =  0. 

A  Mohr-Coulomb  type  yield  model  was  used,  i.e. , 

Y  =  .0003  +  (l-e“P/* 0003)(. 00094  +  1.33  P)  (3) 

where  Y  is  the  yield  strength,  in  megabars.  The  maximum  value  permitted 
far  Y  was  10  kilobars. 
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3.3  CASE  1  -  INTERACTION  OF  STRESS  WA^E  WITH  SINGUB, 

SEMI-INFINITE  CRACK 

The  computational  grid  set  up  for  the  oode  solution  of  Case  1  is 
shown  In  Figure  9.  This  grid  contains  2690  cells,  with  the  basic  cell  size 
set  at  10  cm  x  10  cm.  Beyond  a  central  region  of  interest  the  cell  dimen¬ 
sions  geometrically  increase  in  order  to  conserve  the  total  ceil  count  and 
computational  time.  Representative  results  of  the  code  solution,  as  depicted 
by  particle  velocity  fields  and/or  principal  stress  fields  for  times  of  .3,  .5, 
and  .92  msec,  are  shown  in  Figures  10  to  12  and  in  Figures  2  and  3  in  the 
Summary,  Section  2.3.1.  For  clarity  in  reading  these  plots,  the  field  of 
view  waL  limited  to  the  central  region  of  interest. 

For  the  stress  field  plots,  the  principal  components  of  the  stress 
tensor  for  each  cell  are  shown,  as  follows:  The  magnitude  of  the  two  princi¬ 
pal  stresses  in  the  x-y  plane  are  plotted  in  their  corresponding  principal 
directions.  The  third  principal  stress  (in  the  z  direction)  is  plotted  along 
the  line  bisecting  the  other  two  principal  directions.  Vectors  pointing  to 
the  right  are  compressive,  to  the  left,  tensile.  An  example  of  how  a  stress 
tensor  is  plotted  is  sketched  below: 
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Figure  9.  Initial  Configuration  of  the  Lagrangian  Computational  Grid, 
Case  1 
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Figure  JO.  Particle  Velocity  Field,  Case  1 
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Figure  11.  Particle  Velocity  Field,  Case  1 
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Figure  12.  Principal  Stress  Field,  Case  1 


The  edit*  of  the  velocity  vector  field  plot  the  direction  and  magni¬ 
tude  of  the  velocity  of  each  lattice  point  in  the  computing  grid.  The  vector 
lengths  for  both  the  principal  stress  and  velocity  fields  are  scaled  to  the 
unit  length  Indicated  above  each  plot;  the  units  are  Mb/cm  and  (cm/|isec)/cm , 
respectively.  On  the  more  recently  produced  plots,  the  scale  is  also  graphi¬ 
cally  depicted  in  the  upper  right  hand  corner. 

As  the  stress  wave  interacts  with  the  crack,  the  material  on  the 
right  side  of  the  crack  is  driven  by  the  stress  component  normal  to  the  crack, 
since  shear  stresses  cannot  be  supported.  As  shown  in  Figure  10,  the  velo¬ 
city  vectors  along  the  crack  at  the  shock  front  are  thus  directed  normal  to 
the  crack,  turned  downward  30°.  Also,  as  shown  in  Figure  2,  note  that  the 
principal  stress  tensors  in  the  material  along  the  right  side  of  the  crack  are 
rotated  into  a  direction  transverse  to  the  crack  surface,  again  reflecting  the 
fact  that  the  crack  surface  can  not  bear  shear  stress.  As  the  incident  wave 
runs  along  the  crack,  a  dilatatlonal  wave  and  trailing  shear  wave  are  formed 
which  propagate  across  the  block.  The  transmitted  shock  front  remains 
approximately  planar  and  oriented  at  90°  to  the  x  axis.  In  the  shear  region 
behind  the  shock,  a  distinctly  downward  velocity  flow  is  evident.  This 
action  Induces  material  slippage  along  the  crack,  the  material  on  the  right 
side  of  the  crack  moving  downward  and  to  the  left,  along  the  crack,  relative 
to  the  material  on  the  left  side.  In  addition,  there  was  a  slight  separation, 
or  opening-up,  of  the  materials  on  either  side  of  the  crack,  which  were 
initially  in  contact. 

Time  histories  of  the  material  displacement  at  the  points  indicated 
in  Figure  13  were  recorded  during  the  code  solution.  The  slippage  of 
material  initially  at  the  point  x  =  100  cm,  y  »  0  cm,  as  given  by  the  distance 
between  points  on  opposite  sides  of  the  crack  (points  B  and  C  in  Figure  13), 
is  shown  in  Figure  14,  The  extents  of  the  downward  (y-dlrection)  displace¬ 
ments  of  these  points  are  shown  in  the  time  histories  given  in  Figure  15. 

The  downward  displacements  of  other  points  in  the  field,  at  x  •  55,  158, 

203,  and  253  cm  (points  A,  D,  E,  and  F)  along  the  central  horizontal  plane 
(y  m  0),  are  shown  in  Figure  16.  The  forward  (x-directlon)  displacements  of 
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Figure  16,  Vertical  Displacement  at  Four  Stations  Along  the 
y  *  0  Plane.  Case  1. 
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the«e  points  during  this  time  were  all  ~  0.25  cm.'  The  perturbation  of  flow, 
as  measured  by  the  downward  thrust  of  material,  is  seen  tb  diminish  as  the 
distance  from  the  crack  increases. 

*  ,  i 

Stress  (c^  -  time  profiles  at  points  A,  D,  E,  arid  F  are  shown  in 
Figure  17.  The  peak  stress  in  the  transmitted  wave  (points  D,  E j  and  F)  is. 
seen  to  be  reduced  by  about  25%  from  that  in  the  incident  wave.  The  afore¬ 
mentioned  two-wave  structure  in  the  transmitted  wave  and  the  reflected  wave 
at  point  A  are  also  displayed  in  these  plots.  1 

! 

3.4  CASE  2  -  INTERACTION  OF  STRESS  WAVE  WITH  SINGI£, ' 

FINITE-LENGTH  CRACK  ■  ,  > 

The  initial  configuration  of  the  Lagranglan  grid  set  up  for  this 
problem  is  shown  in  Figure  18.  The  crack  extends  from  the  loading  surface 
(lower  left)  to  the  crack  Up  at  x  ■  200  cm,  y  =  0. 

i 

Representative  results  of  the  SHE P  code  solution  of  this  problem, 
as  depicted  by  the  particle  velocity  fields  for  times  of  .3,  .4,  .5,  ,.6,  and 
1  msec,  are  shown  in  Figures  19  to  22  and  in  Figure  4  ^n  the  Summary,.  1 
Section  2.3.1.  Associated  principal  stress  fields  for  times  of  .4,  .5,  and 
.7  msec  are  shown  in  Figures  23  to  ?5.  Fbr  clarity  in  reading  these  plots, 
the  field  of  view  is  limited  to  the  central  region  of  interest.  1  1 

The  response  of  the  granite  in  this  problem  l!s,  of  necessity, 
similar  to  that  in  the  semi -infinite  crack  problem,  until  the  wave  front 
reaches  the  crack  tip.  Subsequently,  the  wave  system  is  divided  approxi¬ 
mately  in  half,  the  part  above  the  crack  tip  appearing  as  a  simple  plane  1 
wave,  and  that  below  as  a  dilatatiorial  wave  and  trailing  shear  wave,  as  in 
the  previous  solution.  Starting  from  the  cradk  tip,. a  disturbance,  or  bow 
wave,  propagates  into  the  plane  wave  region  above  and  the  "cracked"  * 
region  below,  altering  both  flow  fields  and  creating  an  expanding  region  of 
transition  between  them.  The  diversion  of  flow  around  the  crack  tip  may  be 
seen  in  Figures  19  and  20.  1 
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Figure  18.  Initial  Configuration  of  the  Lagranglan  Computational  Grid, 
Case  2 
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Figure  20.  Particle  Velocity  Field,  Case  2. 


Figure  21.  Particle  Velocity  Field,  Case  2 


Figure  22.  Particle  Velocity  Field,  Case  2. 


Figure  24.  Principal  Stress  Field,  Case  2 
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Figure 2j5.  Principal  Stress  Field,  Case  2. 


Time  histories  of  the  material  displacement  at  the  locations  indi¬ 
cated  in  Figure  26  were  recorded  during  the  code  solution.  The  slippage  of 
material  at  tlxee  points  along  the  crack,  as  given  by  the  distance  between 
points  on  opposite  sides  of  the  crack,  designated  as  (W,X),  (R,S),  and 
(M,N)  In  Figure  26,  is  shown  in  Figure  27.  The  material  on  the  right  side 
of  the  crack  is  moving  downward  and  to  the  left,  along  the  crack,  relative 
to  the  material  on  the  left  side.  The  extents  of  the  downward  (y-dlrectlon) 
displacements  of  these  points  and  the  crack  tip  (point  1)  as  a  function  of 
time  are  shown  in  Figure  28*  The  spatial  trajectory  of  the  point  pair  (W,X) 
on  the  crack  surface  during  the  time  span  of  the  solution  is  plotted  in  Fig¬ 
ure  29.  Time -histories  of  the  vertical  displacements  of  points  above  and 
below  the  crack  tip  along  four  vertical  cuts  (x  m  constant)  through  the  target 
are  shown  in  Figures  30  to  33.  The  downward  shift  of  material  persists  at 
all  these  stations,  but  in  smaller  amounts  as  the  distance  from  the  crack 
Increases. 


Stress  (o^  -  time  profiles  at  points  H,  I,  J,  and  K,  along  the 
horizontal  plane  through  the  crack  tip  (y*  0),  are  shown  in  Figure  34.  Note 
the  Increase  in  stress  over  that  of  the  loading  level  near  the  crack  tip. 

Stress  profiles  along  the  vertical  plane  through  the  crack  tip  (x  m  200  cm) 
are  shown  in  Figure  35.  The  points  above  the  crack  tip  show  a  peak  stress 
of  5  kb  -  corresponding  to  the  loading  or  incident  shock  level,  and  those 
below  the  Up  to  about  4  kb  -  the  transmitted  stress  level.  Stress-Ume  pro¬ 
files  along  a  vertical  plane  to  the  right  of  the  crack  Up,  at  x  m  300  cm,  are 
shown  in  Figure  36.  Here  the  stress  level  is  reduced  at  points  above  as 
well  as  below  the  crack  Up.  Time  histories  of  the  shear  stress  (oxy)  at 
points  along  the  verUcal  plane  through  the  crack  Up  are  shown  in  Figure  37. 

3.5  CASE  3  -  INTERACTION  OF  STRESS  WAVE  WITH  SINGLE, 

FINITE-LENGTH  CRACK,  WITH  CRACK  GROWTH 

For  this  case,  the  same  proolem  as  in  Case  2  was  solved,  but 
with  the  provision  in  the  code  for  permltUng  growth  of  the  crack  acUvated, 
as  described  previously. 
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Figure  39.  Principal  Stress  Field,  Case  3. 


Figure  40.  Principal  Stress  Field,  Case  3. 


Representative  results  of  this  code  solution,  as  shown  by  the 
particle  velocity  fields  for  times  of  .5,  .6,  .7,  and  .9  msec,  are  given  in 
Figures  5  to  7  In  the  Summary,  Section  2.3.1,  and  In  Figure  38.  Associated 
principal  stress  fields  for  times  of  .5  and  .6  msec  are  shown  In  Figures  39 
and  40.  Propagation  of  the  crack  does  occur  for  this  loading;  the  extent  of 
the  crack  at  any  time  Is  Indicated  by  the  circled  lattice  points.  For  this 
loading  function,  the  fracture  criterion  used  was  met  at  relatively  late  times 
after  the  shock  front  has  passed,  so  that  the  effect  of  crack  growth  on  the 
transmitted  wave  Is  not  large. 

4.  ANTI  PLANE  SHEAR  LOADING  OF  A  STATIONARY  CRACK 

4.1  FORMULATION 

The  antiplane  shear  loading  of  a  crack  results  In  equal  and  oppo¬ 
site  out-of-plane  displacements  In  the  quarter  planes  adjacent  to  each  of 
the  two  sides  of  the  crack.  When  this  type  of  shear  loading  Is  applied  to 
a  crack  surface,  no  ln-plane  displacement  components  are  generated.  The 
one  out-of-plane  displacement  component  is  then  governed  by  D*  Alembert's 
equation: 


□w  e  w  -  Aw  *  0  (4) 

s  s 

where 

8  *  ct  (5) 

c  »  (6) 

P 

An  Immediate  advantage  thus  accrues  In  dealing  with  an  out-cf-plane  prob¬ 
lem;  namely  the  solution  can  be  expressed  In  harmonic  functions  rather  than 
blharmonlc  functions.  In  addition,  the  only  elasto-dynanuc  solutions  that 
are  currently  available  for  an  accelerating  crack  are  those  of  Kostrov7, 
Eshelby,8  and  Achenbach9 .  Because  of  these  advantages,  it  was  decided 
to  modify  the  SHEP  code  to  Include  computations  for  the  two  out-of-plane 
shear  stress  components  and  the  one  out-of -plane  displacement  component. 
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fill  three  of  which  depend  only  on  the  in-plane  coordinates, 
equations  are  as  follows: 

The  governing 

do  do  o 

if  I  "  P0® 

(7) 

E  ■  -(P+  q)  V  +  V[(t^+  q^  :  vq  ] 

(8) 

e 

§  (v  q  +  q7) 

(9) 

qd*  f  Q  +  5’) 

(10) 

where 

and 


G 

V 


Is  the  shear  modulus 

Is  the  linear  flctlve  viscosity,  taken  to  be 
proportional  to  the  mesh  size. 


The  difference  form  of  these  equations  follows  the  same  format  as  used  by 
the  rest  of  SHEP.  The  form  of  the  Jauman-Oldroyd  derivative  is  given  in  the 
Appendix. 


4.2  ANALYSIS 

In  order  to  check  out  the  validity  of  the  computer  code,  it  is 
mandatory  to  compare  the  computational  results  with  one  or  more  analytical 
solutions.  The  most  convenient  solution  available  for  this  purpose  is  that 
which  describes  the  propagation  of  an  out-of-plane  (antiplane)  shear  wave 
emanating  from  a  uniformly  loaded  crack.  In  the  following,  the  out-of-plane 
solution  is  presented  and  used  to  check  out  the  corresponding  code  solution 
that  was  performed  for  this  problem. 

Consider  an  infinite  elastic  space  perforated  by  a  semi-infinite 
crack  surface  lying  at  [x  »  0,  -®  <  y  s  0],  To  the  crack  surfaces  x  =  0+and 
x  ■  0  ,  we  apply  a  step  load  in  time  uniformly  distributed,  namely 
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The  plane  [x  =  0,  0  *  y  <  •]  is  a  plane  of  antisymmetry,  on  which  the  out-of- 
plane  displacement  w  is  zero.  In  the  right  half-space  the  displacement  is 
positive,  or  out -of -the -pa  per,  and  vice  versa. 

i 

I  •  : 

This  problem  is  a  standard  one  in  wave  propagation  theory  and  is 
readily  solved  by  the  Wiener-Hopf  technique  and  the  Caigniard  method.  It 
can  also  be  treated  by  the  method  of  characteristics  and, the  theory  of 
Abelian  integrals.  In  either  event,  one  arrives  at  the  solution; 
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l 


»  :  ,1 

wher$  ±  denotes  0i  positive  or  negative  respectively,  superscript  B  denotes 
the  bow  wave  contribution,  and  r  is*  the  radial  distance  from  the  crack  Up. 
The  total  soluUon  is  given  by 


Total 

Bow  Wave 

rf 

Plane  Wave 

G  T 

8  w- 

G  B 
s:  w_  < 

+  (s-x)  H  (s-x)  (13) 

£  WT 

S  + 

i 

i  t 

G  WB  s 

S  + 

(14) 

1 

i  .  j 

4.3  NUMERICAL  SOLUTION 

,  t 

/ 

A  SHEP  code  run  was  set  up  for  this  problem  using  a  grid  containing 
5000  cells  of  uniform  size.1  The  sample  body  was  taken  to  be  1  cm  in  the 
x  dlrecUpn  and'i  1  cm'in  the  y  direcUon,  as  shown  in  the  following  sketch: 
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Since  the  problem  is  antisymmetric,  it  is  sufficient  to  compute  the  wave 
propagation  in  just  one-half  of  the  body  and  Impose  the  boundary  condition 
that  the  displacement  is  zero  along  the  plane  x  =  0,  y  ^  0,  above  the  crack 
tip. 

Q 

The  material  was  characterized  by  a  density  =  2  gm/cm  and  a 
shear  modulus  =  100  kb.  The  crack  was  loaded  uniformly  with  a  shear  stress, 

T)a'  °f  -1  kb- 

With  the  analytic  solution  in  hand,  it  is  convenient  to  perform  a 
number  of  checks  on  the  computational  solution.  We  have  the  following 
derived  relations  for  the  displacement  or  stress  at  three  angles: 


In  addition,  we  have  the  relations: 
PB 


S2c  s 


Pp  -  V  <do  -  s> 


(18) 

(19) 


where  P  denotes  the  rate  of  work  (power  input)  done  on  the  body.  In  Eqns. 
(15  -  17)  we  note  that  the  right-hand  sides  are  all  functions  of  (p).  This 
is  so  because  there  is  no  length-scale  in  the  problem. 
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Figures  41,  42,  and  43  show  respectively  plots  of  Eqn.  (15),  (16), 
and  (17)  versus  —  ,  the  fractional  distance  to  the  wavefront,  along  with  the 
corresponding  results  from  the  code  solution.  Excellent  agreement  is  noted 
in  all  three  cases. 

The  computer  program  also  tabulates  the  total  strain  plus  kinetic 
energy,  the  time  rate  of  change  of  which  is  equal  to  the  power  input.  At  the 
time  the  bow  wave  has  filled  the  sample  body,  the  tabulated  rate  of  energy 
change  was  1.112  x  106  (ergs/cm)/|i  sec,  in  excellent  agreement  with  the 
theoretical  computed  power  input  of  1. 113  x  106  (ergs/cm) /nsec. 

In  Figure  44,  we  show  the  stress  field  at  the  time  at  which  the 
bow  wave  has  just  about  filled  the  sample  body.  At  each  field  point  is 
plotted  a  vector  equal  in  magnitude  to 

M  =  j  <ir>  +  (“s2”)  (20) 

and  having  a  phase  angle  equal  to 

T  vz 

<P  =  arctan  #21) 

yz 

Thus  the  field  vectors  are  given  by: 


V  =  Me1^ 


(22) 


One  can  clearly  see  the  position  of  the  plane  and  bow  waves.  In  the  plot, 
the  front  of  the  plane  wave  appears  to  have  advanced  a  distance  beyond  the 
theoretical  value  of  1.  This  is  due  to  the  non-dissipative  artificial  viscos¬ 
ity  terms  used  in  the  program  which  smoothly  spread  the  wave  across  several 
cells  at  discontinuities,  as  seen,  for  example  in  the  stress  profile  of 
Figure  43.  In  the  bow  wave,  as  the  angle  increases,  the  spreading  is  seen 
to  narrow,  since  the  amplitude  of  the  wave  front  decreases  to  zero  at 
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Figure  41.  Crack-Opening  Displacement  Rrofile 
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Figure  42.  Shear  Stress  Profile  Along  Plane  of  Anti-symmetry  (0  =  3L  ) 
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Figure  44  Shear  Stress  Field. 
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At  the  crack  tip  both  rxz  and  Tyz  become  singular.  The  computer 
run  recorded  highest  values  near  the  crack  tip  (one-half  cell  width  away)  of 
~  5  kb.  The  approach  to  infinity  is  limited  therefore  by  the  fineness  of  the 
mesh  used  in  approximating  the  continuum  by  a  grid.  This,  however,  did 
not  affect  the  good  accuracy  of  the  solution  at  points  in  the  body  away  from 
the  singularity. 

5.  ANTI-PLANE  SHEAR  LOADING  OF  AN  ACCELERATING  CRACK 


5.1  ANALYSIS 


In  the  case  of  a  running  crack,  we  have  an  infinite  elastic  space 
perforated  by  a  semi-infinite  propagation  crack  surface,  lying  at 
[x*  0,  -•  <ys  d  (s)  ],  where  d(s)  is  the  distance  the  crack  has  propaga¬ 
ted  from  the  origin  in  time  s.  The  driving  force  for  the  crack  propagation  is 
provided  by  an  anti -plane  shear  stress 

Txz  ■  -SH(s)  (23) 


applied  to  both  faces  of  the  crack,  and  following  the  crack  as  it  propagates. 
Ahead  of  the  crack,  on  the  plane  x  =  0,  the  displacement  w  is  always 
zero. 


The  above  problem  may  be  attacked  by  first  determining  a  Green's 
displacement  function,  which  is  the  displacement  arising  from  a  unit  force/ 
unit  thickness  applied  at  a  point  on  the  (x  =  0)  surface  designated  by  yQ, 
and  at  a  time  sQ,  for  a  duration  AsQ.  The  solution  to  this  problem  is  well 
known  and  has  the  form 


w 


G 


r*B0 

irG 


H  [s  -  sQ -y/(y-y0)2  +  x2  ] 

•A*  -  sQ)2  -  (y-v0)2  _  x2 


124) 


The  total  displacement  at  any  point  (x,y)  and  at  any  time  s  is  then  given 
by: 
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where 


(25) 

Tq  is  the  stress  applied  ahead  of  the  crack  and  is 
to  be  determined  by  the  condition  that  w  is 
zero  ahead  of  the  crack  on  the  surface  x  =  0. 


The  limits  in  Eqn.  (25)  are  in  part  determined  by  the  Heaviside  function 
appearing  in  Eqn.  (24),  but  are  more  simply  described  after  carrying  out  an 
automorphic  transformation9  to  new  variables. 


* 


_  s_-  y. 

jr 


(26) 


1  = 


s  +  y 


(27) 


These  variables  have  the  property  that  they  allow  the  square  root  appearing 
in  Eqn.  (25)  to  be  factored.  One  can  then  arrive  at  an  Abel  integral  equation 
for  (TQ/S),  under  the  condition  that  w  =  0  ahead  of  the  crack.  After  a  little 
algebra,  we  arrive  at  the  following  results: 


T 

S 


=  —  [arctan 

V 


/s -  y  +  d  (s-y  +  d)  / s  -  y  +  d  (s  -  v  +  d) 

*  y  -  d  ( s  -  y  +  d)  -  *  y  -  d  (s  -  y  +  d) 


] 


(28) 


GD 

S 


where 


—  ^arcsln,/^  (s  +  y.~  'L  +  ^/[d  (s  +  y  -  d)  -  y] [s  +y - d  (s  +y -  d) ] 


(29) 


D  =  wl  _ 


x  =  0 


(30) 


The  function  (d(s-y+d(s-y+d(s-y+.  .  .  has  the  property  that 
d  -•  0  as  y-*  s  and  the  function  d(s+y-d(s+y-d(s+y-.  .  .has  the 
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property  that  d  -*  0  as  y  -»  -s.  We  now  apply  the  Griffith  criterion  (discussed 
below)  /  perform  a  few  operations  with  delta  functions®  ,  and  arrive  at  the 
following  result  for  the  crack  propagation  function: 

~  =(7*  -  1  -  2  *  2  arctan  )  H  (s  -  s*)  (31) 

where 

s*  =  2£G  (32) 

4 

and  where  r  Is  the  specific  surface  energy.  From  Eqn.  (31)  we  conclude 
that 


d 

c 


1-<1T  )2 


(33) 


so  that  the  crack  accelerates  to  sonic  velocity.  The  start  of  crack  propaga¬ 
tion  Is  delayed  by  a  time  s*  until  the  energy  Inside  the  bow  wave  has  built 
up  to  a  critical  value.  Figure  45  shows  how  d/s*  depends  on  s/s* ,  a  no¬ 
parameter  plot. 


5.2  THE  GRIFFITH  CRITERION 

To  improve  the  physical  significance  of  the  numerical  solutions, 
the  conclusion  was  reached  that  a  more  realistic,  dynamic  criterion  for  pre¬ 
dicting  and  following  the  course  of  crack  propagation  in  flawed,  jointed, 
brittle  media  should  be  incorporated  in  the  numerical  method.  A  necessary 
step  in  making  such  a  formulation  change  is  to  check  the  new  code  through 
comparisons  1th  analytical  results  from  a  model  problem,  such  as  the  one 
described  in  the  previous  section. 

The  modifications  needed  to  follow  crack  propagation  are  three¬ 
fold.  These  are  tne  logic  associated  with  the  crack  geometry  and  crack 
prrpagation,  the  velocity  threshold  criterion,  chosen  here  to  be  the  Griffith 
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Figure  45.  Antiplane  Crack  Propagation. 
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criterion,  and  the  velocity  orientation  criterion,  chosen  to  be  the  direction 
transverse  to  the  maximum  principal  stress  vector. 


The  Griffith  criterion  is  a  power  balance.  It  states  that  the  total 
work  done  on  the  body  by  surface  traction  goes  into  increase  of  strain 
energy,  increase  of  kinetic  energy,  and  increase  of  surface  energy.  This 
takes  the  form: 


J*1  P  (f)  d f  =  Ux  +  U2  +  f  (d  -  dQ)  (34) 

o  L 


where 


t* 

P  = 

Ul“ 
U2  = 
W  = 


y  = 


*2  = 


is  a  dummy  time  variable 
<j>  tj  daA  (the  power  input) 

Jw  d  V  (strain  energy) 

J*  Kd  V  (kinetic  energy) 

-2G  Y2  +  \  (K  +  f  G)y2  (strain  energy  density) 
€  £  (dilatation) 

1  O  Ik 

J  (y  -  C k  €  j  )  (second  strain  invariant) 

I  <u);k  +  uk;J>  (straln) 

Ci^  (velocity) 

is  the  stress  tensor 
is  the  displacement  vector 
is  the  tensor  element  of  area 
is  the  volume  element 


H  _  u 

K  =  ~y~  v  v^  (kinetic  energy  density) 

p  is  mass  density 

o 

Uj  and  Ug  are  zero  at  zero  time 
T  is  the  specific  surface  energy 


(35) 

(36) 

(37) 

(38) 

(39) 

(40) 

(41) 

(42) 


(43) 
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d  is  the  Instantaneous  crack  length  (in  2D)  or 
crack  area  (In  3D) 
dQ  Is  the  initial  crack  length 

When  the  crack  Is  not  moving,  d  =0,  in  which  case 

JlP  (f)  df  =  Ux+  U2  (44) 

The  only  basis  for  Eqn.  (44)  becoming  an  inequality  arises  when  the  boun¬ 
dary  conditions  involving  the  crack  surface  change  with  time.  This  state¬ 
ment  requires  elaboration. 

Suppose  we  refer  to  the  problem  discussed  in  the  previous  section, 
namely,  the  one  in  which  shear  stress  S  is  applied  to  part  of  the  boundary 
of  a  rectangle,  the  remaining  boundary  being  clamped.  The  stressed  portion 
of  the  boundary  is  the  crack,  cf.  figure  below 


free 


For  this  case,  when  the  crack  is  not  allowed  to  run,  Eqn.  (44)  is  identic*  lly 
satisfied. 
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i  1 

Imagine  now,  that  after  some  time  t#,  the  crack  starts  to  run,  dQ 
Increases,  and  the  stress  S  is  allowed  to  follow  the  crack.  In  this  case, 
the  left  hand  side  of  Eqn.  (44)  becomes:  i  1 

I 

t*  t.  +At 

f  P(t')  dt'  j  +  f  P(t')  dt'  |  1 

o  d  ®  d  t*  d  =  d+  Ad  1 

o  w  o 


t*  +  2At  1 

+  P  P(t‘)  dt'  |  i  +  .  .  .  .  etc. 

t*+  At  d  =  d  +  2A  d 

w  O  i 


Because  of  the  crack  motion,  energy  is  absorbed  by  the  creation  of  fresh 

•  » 

surface  and  Eqn.  (44)  becomes  an  inequality. 


In  order  to  follow  this  motion,  we  proceed  as  follows:  We  allow 
the  wave  to  propagate  for  some  delay  time  (defined  in  Eqn.  (32)  )  with  the 
crack  at  its  initial  position.  We  then  advance  the  crack  a  small  distance 
Ad,  which  will  be  chosen  to  be  a  fraction  of  cell  length.  We  then  let  the  ' 
wave  propagate  for  one  increment  of  time  At,  and  calculate  the  quantity: 


t*  t  +  A  t 

S  P(f)  df  |  +  J*  P(f)  df  | 

o  d=dn  t*  d  =d  J+Ad 

o  *  o _ 


t  +  At 


-  (Uj  +  u2) 

■  ■  —  —  I 


r 


^calc 


(45) 

t  .  i 

If  this  calculated  Ad  turns  out  to  be  equal  to  the  chosen  Ad*  then  the  time 

I 

increment  At  is  recorded,  and  the  crack  is  advanced  another  Ad.  ,  If  the 
calculated  Ad  turns  out  to  be  less  than  Ad  chosen,  additional  time  incre¬ 
ments  are  allowed  to  elapse  until  equality  is  achieved.  Each1  time  equality 
is  achieved,  the  crack  is  advanced  another  Ad.  A  numerical  plot  of  Ad  vs 
time  can  thu  be  constructed  and  compared  with  the  plot  from  the  analytic 
solution  shown  previously  (Figure  45).  1 
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Ih  setting  up  the  problem  shown  above,  we  have  assumed  that  the 
load  Is  applied  directly  to  the  crack  surface.  In  the  event  a  wave  Is  allowed 
to  impinge  upon  the  crack  by  approach  from  infinity  at  any  angle,  the  problem 
posed  in  tills  fashion  can  always  be  decomposed  into  a  combination  of  a 
plane  Wave  traveling  in  the  absence  of  a  crack,  and  the  problem  shown 
above,  where  the  crbck  is  loaded  directly.  Furthermore  this  problem  can 
always  be  reduced  from  a  full-space  problem  to  a  half-space  problem  by 

I  , 

virtue  of  the  antisymmetry  of  the  shear  loading.  Similarly,  in  dealing  with 
the  in-plane  problem,  we  can  use  symmetry  to  reduce  the  problem  to  that 
of  a  h^lf-space.  Thus  the  problem  posed  above  is  a  very  general  model 
problem.  ! 


Initially,  the  code  was  checked  for, the  case  when  the  crack  is 
not  running  (Eqn.  (44)  ).  To  verify  that  the  computation  of  the  power  integral 
on  tjhe  left-hand  side  of  Eqn.  (44)  matches  within  a  few  percent  of  the  com¬ 
putation  of  ,the  total  energy,  which  is  a  standard  output  of  the  SHEP  program, 

■  i 

it  was  necessary  to  compare  the  computation  of  the  velocity  with  the  ana¬ 
lytical  solution  given  by: 


■  1  ,  =  1  ““  arctn  H  (s-r) 

'  f 

i 

J 

At  the  point  where  the  bow  wave  joins  the  plane  wave,  we  have: 


(46) 


G  w 
Sc 


=  1 


(47) 


Because  of  the  discrete  nature  of  the  computational  algorithm 
used  to  calculate  v&  in  SHEP,  the  computed  velocity  is  not  constant  in 

■  l 

time;  but  rather  fluctuates  around  unity,  the  vibration  gradually  damping  to 
zero  amplitude „  Both  the  amplitude  and  logarithmic  decrement  of  these 
vibrations  can  be  controlled  by  varying  the  artificial  viscosity  and  the  time 
increment.  We  therefore  performed  a  series  of  runs  to  determine  appropriate 
values  for  the  viscosity  coefficient  and  time  step. 
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Figure  46  shows  a  plot  of  normalized  velocity  at  the  wave  front 
as  a  function  of  time  measured  in  numbers  of  cells  in  the  bow  wave.  It  is 
observed  that  critical  damping  occurs  somewhere  between  the  values  of  the 
viscosity  coefficients,  CA  =  .25  and  CA  =  .75.  Variation  of  the  time  step 
had  little  or  no  effect  upon  the  vibration  in  the  range  of  At  studied. 

The  meaning  of  Figure  46  is  that  the  front  of  the  wave  is  spread 
out  over  about  4  cells  on  either  side  of  the  front.  As  we  shall  see  below, 
this  uncertainty  in  &  implies  an  equivalent  uncertainty  in  P,  and,  a  for¬ 
tiori,  an  exaggerated  uncertainty  in  [J**  P(t')  dt'  -  U,  -  U2].  Thus  it  will 

be  extremely  important  to  develop  a  better  form  for  the  fictive  viscosity  in 
order  to  narrow  the  spreading  of  the  wave  front.  Such  a  form  will  be  dis¬ 
cussed  below. 

Holding  in  abeyance  for  the  moment  the  question  of  improving  the 
viscosity  function ,  the  Griffith  criterion  as  displayed  above  was  used  to 
follow  the  propagation  of  the  semi-infinite  crack.  The  wave  was  allowed  to 
propagate  for  90  cycles  with  the  crack  unextended;  then  the  crack  was 
extended  one  cell,  and  the  additional  time  needed  for  dynamic  instability 
to  arise  was  computed  according  to  the  algorithm  provided  above.  In  this 
way  the  crack  was  extended  three  times ,  and  the  results  of  the  computation 
are  shown  in  Figure  47.  We  note  that  the  crack  runs  effectively  slower  than 
the  theoretical  sonic  value.  Further  work  remains  to  be  done  to  refine  the 
algorithm  and  the  viscosity. 

5.3  THE  FICTIVE  VISCOSITY 

As  soon  as  a  continuum  is  replaced  by  a  discrete  grid,  it  becomes 
possible  to  replace  a  slowly  changing  function  by  a  periodic  one.  For 
example,  the  simple  function 

y=  0  (48) 

can  be  replaced  by 
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Figure  47.  Crack  Propagation  History 


y  =  A  sin  nir 


(49) 

where 


This  Is  an  extreme  case  but  it  demonstrates  the  point.  What  the  difference 
equations  actually  offer  depends  on  the  outcome  of  a  stability  analysis. 

Such  analyses  have  been  carried  out  and  they  lead  to  the  idea  that  a  flctive 
viscosity  must  be  incorporated  in  the  constitutive  law  to  produce  damping 
and  stability.  The  penalty  that  one  pays  is  a  spreading  of  a  wave  front. 

Von  Neumann  and  Richtmyer10  suggested  that  a  quadratic  vis¬ 
cosity  term  be  used  to  produce  effective  damping,  and  this  has  indeed  been 
proved  to  be  a  useful  tool.  Other  people  have  used  a  linear  viscosity  term 
with  less  striking  results.  The  linear  viscosity  spreads  the  wave  in  an 
increasing  manner.  Thus  it  appears  that  the  quadratic  viscosity  is  an 
Indispensable  tool. 

For  the  out-of-plane  motion  that  we  have  been  discussing,  there 
is  currently  not  available  a  quadratic  viscosity  formulation;  thus  a  linear 
viscosity  was  used  in  the  computation  discussed  above.  It  Is  possible 
however  to  suggest  a  cubic  formulation,  and  the  formulation  which  we  pro¬ 
pose  below  will  be  evaluated  during  work  on  a  related  program.  The  etymo¬ 
logy  is  as  follows. 

The  constitutive  theory  of  fluid  dynamics  leads  to  the  most  general 
representation  of  one  isotropic  tensor  (the  stress  tensor)  in  terms  of  another 
isotopic  tensor  (the  strain  rate  tensor),  namely; 

t'  a  -  pF  +  ajH  +a2  H  .  d  (51) 

where  p  is  the  hydrostatic  pressure  and  and  a2  are  scaler  functions  of 
the  scaler  invariants;  {  Id,  lld ,  md } 
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where 


(52) 


Id  =  Tr(5)=7^=J  =  "  £ 

Hd  =  IIId  Tx  {  d  ~l)  (53) 

md  =  Det{5)  (54) 

The  terms  multiplied  by  (ttj,  a2)  represent  the  dissipative  contributions  to 
the  total  stress.  One  way  of  incorporating  quadratic  dissipation  is  to  take 
[a2  -  constant].  This  term  however  contributes  only  to  the  normal  stresses, 
not  to  the  shear  stresses.  The  result  produces  what  is  known  as  a  "normal 
stress"  or  Welssenberg  effect.  This  effect  is  not  useful  for  the  flows  that 
we  deal  with  in  fracture  mechanics. 

An  alternative  way  of  introducing  nonlinear  viscosity  is  to  set 
[a2  =  0]  and  take: 

«1=  2voP  (1  -r  Id)  (55) 

where  vQ  is  the  kinematic  viscosity 

and  r  Is  a  relaxation  time. 

This  form  is  very  interesting  for  several  reasons.  First  of  all,  the  product 
(p  vQ)  is  the  linear  viscosity  coefficient  r)Q,  Thus  the  effect  of  linear  vis¬ 
cosity  is  included.  Because  of  conservation  of  mass,  Eqn.  (55)  may  be 
rewritten  as: 


«1  =  2  VQ  (p  +  T  p)  (56) 

which  evinces  the  typical  form  of  viscous  relaxation.  Thirdly  the  quadratic 
term  (rld)  enters  in  the  form  according  to  which  real  materials  behave; 
namely,  the  apparent  viscosity  decreases  with  increasing  shear  rate. 
Finally,  if  the  lineai  term  in  Eqn.  (55)  is  dropped,  the  expression  above 
reduces  exactly  to  the  form  suggested  by  Richtmeyer  and  Von  Neumann. 
Thus  Eqn.  (56)  displays  a  useful  form  for  introducing  quadratic  viscosity 
into  all  components  of  th9  stress  tensor. 
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In  the  event  the  flow  is  isochoric  as  it  is  in  the  case  of  anti- 
plane  shear,  one  can  not  introduce  a  quadratic  viscosity,  because  1^=0. 
One  can  however  introduce  a  cubic  viscosity  by  using  the  form: 

«1=  2  vop  (1+  ;3lld)  (57) 

In  both  forms  suggested  forft^,  vQ  can  be  taken  proportional  to  the  square 
of  the  mesh  size.  In  the  compressible  form,  involving  quadratic  viscosity, 
r  can  be  taken  proportional  to  the  time  step,  and  in  the  isochoric  form, 
involving  cubic  viscosity,  /?  can  be  taken  proportional  to  the  square  of  the 
time  step. 


For  the  compressible  or  quadratic  form,  we  have  repeated  the 
Von  Neumann-Rlchtmeyer  analysis.  This  analysis  leads  to  the  equation: 


dx  =  -2  r  du 

TT+fr*  i)  u  (u_  -  u)  - 1 

vo  p 


(58) 


It  is  easy  to  show  that,  as  u-  {0,  up}  ,  x  -*  .  This  is  in  major  contra¬ 

diction  to  the  Von  Neumann  analysis,  where  x  -»  a  finite  length.  What  Von 
Neumann  and  Richtmeyer  then  did  was  to  tack  on  to  the  shock  front  two 
straight-line  constant  velocity  regions,  with  a  concatenate  discontinuity 
in  curvature.  Because  the  flow  equations  are  second  order,  this  solution 
does  not  satisfy  the  original  differential  equation  and  must  be  regarded  as 
an  approximate  representation.  The  representation  shown  above  does  not 
suffer  from  this  limitation,  and  can,  by  variation  of  r  and  v  ,  be  made  to 
approach  the  shock  discontinuity  as  closely  as  desired.  Figure  48  shows 
a  plot  of  particle  velocity  profile  for 


~  =  10"5,  10" 6,  and  10"7 
vo 

-8 

and  t  =  10  sec. 

One  observes  that  the  shock  thickness  can  be  made  to  take  on  a  broad  range 
of  values. 
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Figure  48.  Particle  Velocity  Profile 
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APPENDIX 


CORRECTION  TO  THE  ROTATIONAL  TERM  IN  THE  STRESS  CALCULATIONS 

During  the  examination  of  the  formulation  of  the  code  preparatory  to 
making  the  modifications  discussed  in  Sections  4  and  5  of  the  report,  an  error 
in  the  sign  of  the  rotational  correction  term  in  the  stress  calculations,  as 
originally  proposed  by  Wilkins^*,  was  discovered.  In  addition,  from  the 

hO 

standpoint  of  material  objectivity  ,  the  form  of  the  rotational  term  is  seen 
to  be  written  in  only  an  approximate  form.  These  statements  are  documented 
below.  The  sign  was  corrected  in  our  code  and  a  check  of  the  correction 
was  made  by  running  a  model  problem.  The  results  confirmed  the  sign  correc¬ 
tion  and,  in  addition,  verified  the  accuracy  of  the  SHEP  code. 

In  SHEP,  the  constitutive  equations  are  correctly*  subsumed  in  the 

form  : 


■  su"6ij 

p  0) 

(Al) 

% 

■  2g  Kj 

-i  *y 

(A2) 

e 

=  InJ 

(A3) 

J 

dV 

(A4) 

where 

tjj  is  the  cartesian  Couchy  stress  tensor 

s^  is  the  stress  de  via  tor  tensor 
P  is  the  pressure 

J  is  the  local  volume  ratio,  equal  to  the  square 
root  of  tha  third  stretch  invariant. 

♦The  term  "correctly"  implies  that  Eqn.  (A2)  satisfies  the  principle  of 
material  objectivity. 
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is  the  shear  modulus 


6 

6  is  the  dilatation 

is  the  Cartesian  strain  rate  tensor 


and  in  (A4) , 


is  the  local  instantaneous  specific  volume 
is  the  particle-differentiator. 


From  Eqn.  (A3)  and  (A4),  it  follows  that 

J  =  Y  (A5) 

Furthermore  in  Eqn.  (A2),  the  dot  over  e^  denotes  the  material  derivative 
while  the  hat  over  s^  denotes  the  Jaumann-Oldroyd  derivative.  The  form 
of  the  latter  is  given  by: 


where 


or 


where 


hk  "  "ik  -  vi,J  *jk  +  *i)  vJ,k 
V!,J  *  +  “lj 


(Ojj  is  the  spin  tensor. 


(A6) 

(A7) 

(A8) 


The  last  two  terms  on  the  right  hand  side  of  Eqn.  (A6)  correspond  to 
Wilkins'  6jjj  a  simple  calculation  reveals  that  the  entire  form  of  Eqn.  (B-14) 
on  page  78  of  Reference  A1  is  formally  incorrect  because  of  a  sign  error. 


The  sign  error  apparently  arose  in  the  following  way:  Wilkins  intro¬ 
duces  an  angular  velocity  vector,  which  he  denotes  by  [sin  oo  ].  The  angular 
velocity  vector  is  correctly  taken  as  the  curl  of  the  velocity  vector  which, 
according  to  the  usual  right-handed  correction,  corresponds  to  a  counter- 
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clockwise  rotation.  Wilkins  then  uses  the  rotation  angle  to  to  transform  the 


stress  components  to  new  rotated  components.  The  transformation  equations 

A3 

he  adopts  are  taken  from  Timoshenko  .  Presumably  Wilkins  changed  the 
signs  of  the  terms  that  are  odd  in  to,  because  Timoshenko  shows  a  diagram 
in  which  the  rotation  is  clockwise.  What  one  must  notice,  however,  is  that 
Timoshenko's  picture  is  based  on  a  left-handed  system  of  axes.  Thus  Timo¬ 
shenko's  equation  should  have  been  subsumed  without  the  sign  change. 

To  further  check  our  logic,  we  ran  the  following  simple  problem.  A 
slab  is  set  into  simple  shear  motion,  cf: 


Nonlinear  (Poynting)  effects  are  associated  with  the  ,’engthening  of  the 
originally  upright  fibers.  The  correct  rotation  terms  (with  the  right  sign)  pre¬ 
dict  a  tensile  stress  in  the  lengthened  fibers,  in  agreement  with  the  analytical 
solution.  The  original  program,  with  the  incorrect  sign,  predicts  compression. 
Keep  in  mind  that  these  are  second  order  stresses  and  thus  the  effect  of  this 
error  on  the  overall  computation  becomes  important  only  when  the  fhst  order 
stresses  are  comparable  to  the  shear  modulus  of  the  material.  The  reason  is 
that  the  second-order  stresses  go  as  the  square  of  the  first-order  stresses, 
e.g. ,  in  the  shear  problem: 

1  ~  >2  <»> 

\ 

In  the  process  of  making  these  chang is  in  SHEP,  we  checked  out  a 
model  problem,  namely  the  simple  shear  at  constant  velocities  of  an  infinitely 
long  slab  (a  1-D  problem).  The  results,  in  terms  of  the  stress  at  the  moving 
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surface  and  the  stress  at  the  fixed  surface  are  plotted  in  Figures  A1  and  A2, 
respectively.  The  analytical  solution  is  just  a  sum  of  Heaviside  step  functions. 
The  agreement  with  the  analytical  solution  is  attested  to  by  the  fact  that  the 
computed  jumps  fall  directly  (within  1%)  on  the  grid  line  (both  coordinate  axes 
are  non-dimensionalized).  The  general  solution  looks  like 


where 


T  » 

"  IjFg  =  S  [H(s-y-2nh]+  H  (s  +  y  -  2h  -  2nh)  ] 

v 

MQ  (the  Mach  number)  = 
s  =  ct 


and  h  is  the  height  of  the  slab. 


(A10) 


I 
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Figure  Al.  Shear  Stress  at  the  Moving  Surface 
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Figure  A2.  Shear  Stress  at  the  Clamped  Surface 
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